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Neural spike trains, which are sequences of very brief jumps in 

voltage across the cell membrane, were one of the motivating ap- 

^^ ■ plications for the development of point process methodology. Early 

work required the assumption of stationarity, but contemporary ex- 
periments often use time-varying stimuli and produce time-varying 
neural responses. More recently, many statistical methods have been 
'^^ \ developed for nonstationary neural point process data. There has also 

been much interest in identifying synchrony, meaning events across 
two or more neurons that are nearly simultaneous at the time scale of 
C/3 . the recordings. A natural statistical approach is to discretize time, us- 

ing short time bins, and to introduce loglinear models for dependency 
among neurons, but previous use of loglinear modeling technology 
^ ' has assumed stationarity. We introduce a succinct yet powerful class 

^N) I of time-varying loglinear models by (a) allowing individual-neuron 

C""^ ■ effects (main effects) to involve time-varying intensities; (b) also al- 

OO ' lowing the individual-neuron effects to involve autocovariation effects 

(history effects) due to past spiking, (c) assuming excess synchrony 
f^^ ' effects (interaction effects) do not depend on history, and (d) as- 

^^ , suming all effects vary smoothly across time. Using data from the 

primary visual cortex of an anesthetized monkey, we give two exam- 
ples in which the rate of synchronous spiking cannot be explained 
by stimulus-related changes in individual-neuron effects. In one ex- 
ij ■ ample, the excess synchrony disappears when slow- wave "up" states 

I ^N , are taken into account as history effects, while in the second example 

j^ ' it does not. Standard point process theory explicitly rules out syn- 

chronous events. To justify our use of continuous-time methodology, 
we introduce a framework that incorporates synchronous events and 
provides continuous-time loglinear point process approximations to 
discrete-time loglinear models. 
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1. Introduction. One of the most important tecimiques in learning about 
the functioning of the brain has involved examining neuronal activity in lab- 
oratory animals under varying experimental conditions. Neural information 
is represented and communicated through series of action potentials, or spike 
trains, and the central scientific issue in many studies concerns the physi- 
ological significance that should be attached to a particular neuron firing 
pattern in a particular part of the brain. In addition, a major relatively new 
effort in neurophysiology involves the use of multielectrode recording, in 
which responses from dozens of neurons are recorded simultaneously. Much 
current research focuses on the information that may be contained in the 
interactions among neurons. Of particular interest are spiking events that 
occur across neurons in close temporal proximity, within or near the typical 
one millisecond accuracy of the recording devices. In this paper we provide 
a point process framework for analyzing such nearly synchronous events. 

The use of point processes to describe and analyze spike train data has 
been one of the major contributions of statistics to neuroscience. On the one 
hand, the observation that individual point processes may be considered, ap- 
proximately, to be binary time series allows methods associated with gener- 
alized linear models to be applied [cf. Brillinger (1988, 1992)]. On the other 
hand, basic point process methodology coming from the continuous-time 
representation is important both conceptually and in deriving data-analytic 
techniques [e.g., the time-rescaling theorem may be used for goodness of fit 
and efficient spike train simulation; see Brown et al. (2001)]. The ability to go 
back and forth between continuous time, where neuroscience and statistical 
theory reside, and discrete time, where measurements are made and data are 
analyzed, is central to statistical analysis of spike trains. From the discrete- 
time perspective, when multiple spike trains are considered simultaneously 
it becomes natural to introduce loglinear models [cf. Martignon et al. (2000)] 
and a widely read and hotly debated report by Schneidman et al. (2006) ex- 
amined the extent to which pairwise dependence among neurons can capture 
stimulus-related information. A fundamental limitation of much of the work 
in this direction, however, is its reliance on stationarity. The main purpose 
of the framework described below is to handle the nonstationarity inherent 
in stimulus-response experiments by introducing appropriate loglinear mod- 
els while also allowing passage to a continuous-time limit. The methods laid 
out here are in the spirit of Ventura, Cai and Kass (2005), who proposed 
a bootstrap test of time-varying synchrony, but our methods are different in 
detail and our framework is much more general. 

Statistical modeling of point process data focuses on intensity functions, 
which represent the rate at which the events occur, and often involve covari- 
ates [cf. Brown et al. (2004), Kass, Ventura and Brown (2005), Paninski et 
al. (2009) and references therein]. A basic distinction is that of conditional 
versus marginal intensities: the conditional intensity determines the event 
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Fig. 1. Neural spike train raster plots for repeated presentations of a drifting sine wave 
grating stimulus. (A) Single cell responses to 120 repeats of a 10 second movie. At the top 
is a raster corresponding to the spike times, and below is a pen-stimulus time histogram 
(PSTH) for the same data. Portions of the stimulus eliciting firing are apparent. (B) The 
same plots as m (A), for a different cell. (C) Population responses to the same stimulus, 
for 5 repeats. Each block, corresponding to a single trial, is the population raster for v = 128 
units. On each trial there are several dark bands, which constitute bursts of network activity 
sometimes called "up states. " Up state epochs vary across trials, indicating they are not 
locked to the stimulus. 



rate for a given realization of the process, while the marginal intensity is the 
expectation of the conditional intensity across realizations. In neurophysio- 
logical experiments stimuli are often presented repeatedly across many tri- 
als, resulting in many replications of the multiple sequences of spike trains. 
This is the situation we concern ourselves with here, and it is illustrated in 
Figure 1, part A, where the responses of a single neuron for 120 trials are 
displayed: each line of the raster plot shows a single spike train, which is 
the neural response on a single trial. The experiment that generated these 
data is described in Section 1.1. The bottom panel in part A of Figure 1 dis- 
plays a smoothed peristimulus time histogram (PSTH), which summarizes 
the trial- averaged response by pooling across trials. As we explain in greater 
detail in Section 1.2, scientific questions and statistical analyses may con- 
cern either within-trial responses (conditional intensities) or trial-averaged 
responses (marginal intensities). 

A point process evolves in continuous time but, as we have noted, it is 
convenient for many statistical purposes to consider a discretized version. 
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Decomposing time into bins of width 6, we may define a binary time series 
to be 1 for every time bin in which an event occurs, and for every bin 
in which an event does not occur. It is not hard to show that, under the 
usual regularity condition that events occur discretely (i.e., no two events 
occur at the same time), the likelihood function of the binary time series 
approximates the likelihood of the point process as (5 — )■ 0. For a pair of 
point processes, the discretized process is a time series of 2 x 2 polytomous 
variables indicating, in each time bin, whether an event of the first process 
occurred, an event of the second process occurred, or both, or neither. This 
suggests analyzing nearly synchronous events based on a loglinear model 
with cell probabilities that vary across time. Intuitive as such procedures 
may be, their point process justification is subtle: the standard regularity 
condition forbids two processes having synchronous events, so it is not ob- 
vious how we might obtain convergence to a point process (as (5 — > 0) for 
discrete-process likelihoods that incorporate synchrony. 

One way out of this impasse is to introduce a marked point process frame- 
work in which each event/mark could be of three distinct types: first process, 
second process, or both. The standard marked point process requires modi- 
fication, however, because it fails to accommodate independence as a special 
case. Under independence, the discretized events for each process occur with 
probability of order 0{5), while the synchronous events occur with proba- 
bility of order 0((5^) as (5 — t- 0. We refer to this as a sparsity condition, and 
the generalization to multiple processes involves a hierarchical sparsity con- 
dition. Once we introduce a family of marked point processes indexed by 5, 
we can guarantee hierarchical sparsity. Not only does this allow, as it must, 
the special case of independence models, but it also makes the conditional 
intensity for neuron i depend only on the history for neuron i, asymptoti- 
cally (as (5 — 7> 0). This in turn avoids confounding the dependence described 
by the loglinear model and greatly reduces the dimensionality of the prob- 
lem. We require two very natural regularity conditions based on well-known 
neurophysiology: the existence of a refractory period, during which the neu- 
ron cannot spike again, and smoothness of the conditional intensity across 
time. It would be possible, and sometimes advantageous, instead to model 
dependence through the individual- neuron conditional intensity functions. 
The loglinear modeling approach used here avoids this step. 

1.1. A motivating example. In a series of experiments performed by one 
of us (Kelly, together with Dr. Matthew Smith), visual images were dis- 
played at resolution 1024 x 768 pixels on a computer monitor, while the 
neural responses in the primary visual cortex of an anesthetized monkey 
were recorded. Each of 98 distinct images consisted of a sinusoidal grating 
that drifted in a particular direction for 300 milliseconds, and each was re- 
peated 120 times. Each repetition of the complete sequence of stimuli lasted 
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approximately 30 seconds. This kind of stimulus has been known to drive 
cells in the primary visual cortex since the Nobel prize-winning work of 
Hubel and Wiesel in the 1960s. With improved technology and advanced 
analytical strategies, much more precise descriptions of neural response are 
now possible. A small portion of the data from 5 repetitions of many stimuli 
is shown in part C of Figure 1 . 

The details of the experiment and recording technique are reported in 
Kelly et al. (2007). A total of 125 neural "units" were obtained, which in- 
cluded about 60 well-isolated individual neurons; the remainder were of un- 
determined origin (some mix of 1 or more neurons) . The goal was to discover 
the interactions among these units in response to the stimuli. Each neuron 
will have its own consistent pattern of responses to stimuli, as illustrated 
in parts A and B of Figure 1. Synchronous spiking across neurons is rela- 
tively rare. However, in each of the 5 blocks within part C of Figure 1 (each 
block corresponding to a single trial) several dark bands of activity across 
most neurons may be seen during the trial. These bands correspond to what 
are often called network "up" states, and are often seen under anesthesia. 
For discussion and references see Kelly et al. (2010). It would be of inter- 
est to separate the effects of such network activity from other synchronous 
activity, especially stimulus-related synchronous activity. The framework in 
this paper provides a foundation for statistical methods that can solve such 
problems. 

1.2. Overview of approach. We begin with some notation. Suppose we 
observe the activity of an ensemble of v neurons labeled 1 to ly over a time 
interval [0,r), where T > is a constant. Let A'^ denote the total number 
of spikes produced by neuron i on [0,T) where i = l,...,z/. The resulting 
(stochastic) sequence of spike times is written as < Sj^ < • • • < s* ^ < T. For 

the moment we focus on the case z/ = 3, but other values of ;/ are of interest 
and with contemporary recording technology z^ ~ 100 is not uncommon, as 
in the experiment in Section 1.1. Let 5 > be a constant such that T is 
a multiple of 5 (for simplicity). We divide the time interval into bins of 
width 6. Define X^{t) = 1 if neuron i has a spike in the time bin [t,t + 6) 
and otherwise. Because of the existence of a refractory period for each 
neuron, there can be at most 1 spike in [t,t + 6) from the same neuron if 5 
is sufficiently small. Then writing 

P^Mit) = P{X\t) = a,X\t) = b,X\t) = c) ya,b,ce {0, 1}, 

the data would involve spike counts across trials [e.g., the number of trials 
on which (X^(t), X^{t),X^{t)) = (1,1,1)]. The obvious statistical tool for 
analyzing spiking dependence is loglinear modeling and associated method- 
ology 
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Three complications make the problem challenging, at least in principle. 
First, there is nonstationarity: the probabilities vary across time. The data 
thus form a sequence of 2^ contingency tables. Second, absent from the above 
notation is a possible dependence on spiking history. Such dependence is 
the rule rather than the exception. Let 'HJ denote the set of values of X*(s), 
where s < t, and s, t are multiples of b. Thus, l-Lt = {T-L\ , . . . , l-ii) is the history 
of the binned spike train up to time t. We may wish to consider conditional 
probabilities such as 

PliiiA'fit) = P{X\t) = a,X\t) = b, X^{t) = c\nt) 

for a,b,c£ {0,1}. Third, there is the possibility of precisely timed lagged 
dependence (or time-delayed synchrony): for example, we may want to con- 
sider the probability 

(1) Pl;^;^{s,t,u) = p{x\s) = i,x\t) = i,x'{u) = 1), 

where s,t,u may be distinct. Similarly, we might consider the conditional 
probability 

pl;^;^{s,t,u\nl,nlnl) = p{x\s) = i,x\t) = i,x^{u) = i\nl,nlnl). 

In principle, we would want to consider all possible combinations of lags. 
Even for z^ = 3 neurons, but especially as we contemplate z^ ^ 3, strong re- 
strictions must be imposed in order to have any hope of estimating all these 
probabilities from relatively sparse data in a small number of repeated trials. 
To reduce model dimensionality, we suggest four seemingly reasonable tac- 
tics: (i) considering models with only low-order interactions, (ii) assuming 

the probabilities PA ^ (t) or P^ '^^ ^ (^l^t) vary smoothly across time t, (iii) re- 
stricting history effects to those that modify a neuron's spiking behavior 
based on its own past spiking, and then (iv) applying analogues to standard 
loglinear model methodology. Combining these, we obtain tractable models 
for multiple binary time series to which standard methodology, such as max- 
imum likelihood and smoothing, may be applied. In modeling synchronous 
spiking events as loglinear time series, however, it would be highly desir- 
able to have a continuous-time representation, where binning becomes an 
acknowledged approximation. We therefore also provide a theoretical point 
process foundation for the discrete multivariate methods proposed here. 

It is important to distinguish the probabilities -P^'^^ {t) and -P^'^^ (tlT-Lt). 
The former are trial- averaged or marginal probabilities, while the latter 
are within-trial or conditional probabilities. Both might be of interest but 
they quantify different things. As an extreme example suppose, as some- 
times is observed, each of two neurons has highly rhythmic spiking at an 
approximately constant phase relationship with an oscillatory potential pro- 
duced by some large network of cells. Marginally these neurons will show 
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strongly dependent spiking. On the other hand, after taking account of the 
oscillatory rhythm by conditioning on each neuron's spiking history and/or 
a suitable within-trial time-varying covariate, that dependence may vanish. 
Such a finding would be informative, as it would clearly indicate the nature 
of the dependence between the neurons. In Section 3 we give a less dramatic 
but similar example taken from the data described in Section 1.1. 

We treat marginal and conditional analyses separately. Our use of two 
distinct frameworks is a consequence of the way time resolution will affect 
continuous-time approximations. We might begin by imagining the situation 
in which event times could be determined with infinite precision. In this case 
it is natural to assume, as is common in the point process literature, that 
no two processes have simultaneous events. As we indicate, this conception 
may be applied to marginal analysis. However, the event times are necessar- 
ily recorded to fixed accuracy, which becomes the minimal value of S, and 5 
may be sufficiently large that simultaneous events become a practical pos- 
sibility. Many recording devices, for example, store neural spike event times 
with an accuracy of 1 millisecond. Furthermore, the time scale of physiolog- 
ical synchrony — the proximity of spike events thought to be physiologically 
meaningful — is often considered to be on the order of 5 = 5 milliseconds [cf. 
Griin, Diesmann and Aertsen (2002a, 2002b) and Griin (2009)] . For within- 
trial analyses of synchrony, the theoretical conception of simultaneous (or 
synchronous) spikes across multiple trials therefore becomes important and 
leads us to the formalism detailed below. The framework we consider here 
provides one way of capturing the notion that events within 6 milliseconds 
of each other are essentially synchronous. 

The rest of this article is organized as follows. Section 2 presents the 
methodology in three subsections: Sections 2.1 and 2.2 introduce marginal 
and conditional methods in the simplest case, while Section 2.3 discusses the 
use of loglinear models and associated methodology for analyzing spiking 
dependence. In Section 3 we illustrate the methodology by returning to 
the example of Section 1.1. The main purpose of our approach is to allow 
covariates to take account of such things as the irregular network rhythm 
displayed in Figure 1, so that synchrony can be understood as either related 
to the network effects or unrelated. Figure 2 displays synchronous spiking 
events for two different pairs of neurons, together with accompanying fits 
from continuous-time loglinear models. For both pairs the independence 
model fails to account for synchronous spiking. However, for one pair the 
apparent excess synchrony disappears when history effects are included in 
the loglinear model, while in the other pair they do not, leading to the 
conclusion that in the second case the excess synchrony must have some 
other source. Theory is presented in Sections 4-6. We add some discussion 
in Section 7. All proofs in this article are deferred to the Appendix. 
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Fig. 2. Synchronous spike analysis for two pairs of neurons. Results for one pair shown 
on left, in parts (A)-(D) and for the other pair on the right in parts (E)-(H). Part (A) 
Response of a cell to repetitions of a 1 second drifting grating stimulus. The raster plot 
is shown above and the smoothed PSTH below. Part (B) Response from a second cell, as 
in (A). In both (A) and (B), spikes that are synchronous between the pair are circled. 
Part (C) Correct joint spike predictions from model, shown as circles [as m parts (A) 
and (B)y, when false positive rate is set at 10%. In top plot the joint spikes are from 
the history-independent model, as in (13), while in the bottom plot they are as in (16), 
including the network covariate in the history term. Part (D) ROC curves for the models 
m part (C). Parts (E), (F), (G) and (H) are similar to Parts (A), (B), (C) and (D) but 
for the second pair of neurons. 
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2. Methodology. In this section we present our continuous-time loghnear 
modehng methodology. We begin with the simplest case of z/ = 2 neurons, 
presenting the main ideas in Sections 2.1 and 2.2 for the marginal and condi- 
tional cases, respectively, in terms of the probabilities -P^ (i) = P{X^{t) = a), 
Plf{t) = P{X^{t) = a,X'^{t) = b), etc., for aU a,be {0, 1}. We show how we 
wish to pass to the continuous-time limit, thereby introducing point pro- 
cess technology and making sense of continuous-time smoothing, which is 
an essential feature of our approach. In Section 2.3 we reformulate using 
loglinear models, and then give continuous-time loglinear models for z^ = 3. 
Our analyses in Section 3 are confined to i^ = 2 and z^ = 3 because of the 
paucity of higher-order synchronous spikes in our data. Our explicit models 
for z^ = 3 should make clear how higher-order models are created. We give 
general recursive formulas in Sections 5 and 6. 

2.1. Marginal methods for v = 2. The null hypothesis 

(2) H^ : pI'^ {t) = Pi {t)Pl {t) Vt G r = {0, 5, 25, . . . , T - <5} 

is a statement that both neurons spike in the interval [t, t + J), on the average, 
at the rate determined by independence. Defining C,{t) by 

(3) Pl'^{t)=Pl{t)Pl{t)at) yteT, 
we may rewrite (2) as 

(4) Ho:C{t) = i yter. 

As in Ventura, Cai and Kass (2005), to assess Hq, the general strategy we 
follow is to (i) smooth the observed-frequency estimates of P^it), -Pf(t) 
and P^ \ (t) across time t, and then (ii) form a suitable test statistic and 
compute a p-value using a bootstrap procedure. We may deal with time- 
lagged hypotheses similarly, for example, for a lag /i > 0, we write 

Pl'^{t,t + 6h) = P{X^it) = l,X^{t + 6h) = l) 

(5) 

= Pl{t)P^{t + 6h)C{t,t + 6h), 

then smooth the observed-frequency estimates for P^ \ {t,t + 6h) as a function 
of t, form an analogous test statistic and find a p-value. 

To formalize this approach, we consider counting processes N^ correspond- 
ing to the point processes S]^,S2, . . . ,s* i, i = 1,2 (as in Section 1.2 with 

u = 2). Under regularity conditions, the following limits exist: 
(6) Ai'2(t) = hm 5-2p((iV,V - Nl)iN^^, - N?) = 1), 

o— >U 

e(t) = limC(t). 
5— >-0 
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Consequently, for small 6, we have 



Pm^X\t)6, Pl;t{t)^y^'{t)6 



The smoothing of the observed-frequency estimates for P^ \ (t) may be un- 
derstood as producing an estimate A^'^(t) for A^'^(t). The null hypothesis 
in (2) becomes, in the limit as 6—^0, 

Ho:X^''^{t) = X\t)X^{t) Vte[0,T), 

or, equivalently, 

(7) Ho:m = ^ VtG[0,r). 

The lag h case is treated similarly. Under mild conditions. Theorems 1 and 2 
of Section 5 show that the above heuristic arguments hold for a continuous- 
time regular marked point process. This in turn gives a rigorous asymptotic 
justification (as 6—^0) for estimation and testing procedures such as those in 
steps (i) and (ii) mentioned above, following (4), and illustrated in Section 3. 

2.2. Conditional methods for u = 2. To deal with history effects, equa- 
tion (3) is replaced with 

(8) p^mt) = pimDplmDm vt g r, 

where i-Ll, i = 1,2, are, as in Section 1, the binned spiking histories of neu- 
rons 1 and 2, respectively, on the interval [0,i). Analogous to (4), the null 
hypothesis is 

Hq : c(t) = 1 vt G r. 

We note that there are two substantial simplifications in (8). First, P{(t|?^() = 
P^tlTil), which says that neuron i's own history Til is relevant in modifying 
its spiking probability (but not the other neuron's history). Second, (^{t) does 
not depend on the spiking history "Hf This is important for what it claims 
about the physiology, for the way it simplifies statistical analysis, and for 
the constraint it places on the point process framework. Physiologically, it 
decomposes excess spiking into history-related effects and stimulus-related 
effects, which allows the kind of interpretation alluded to in Section 1 and 
presented in our data analysis in Section 3. Statistically, it improves power 
because tests of Hq effectively pool information across trials, thereby in- 
creasing the effective sample size. 

Consider counting processes N^, i = 1, 2, as in Section 2.1. Under regular- 
ity conditions, the following limits exist for t G [0,r): 

X\t\nD = ^m5~^P{Ni^s - N'i = im, i = 1, 2, 
(9) X''\t\nt) = hm r2p((iV,V5 - Nl){N^^s - N?) = M-Ht), 

0— >U 

at) = limC(t), 
5^0 
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where T-Lt = lim^^o^i and T-Ll = lim^^o^t)^ = 1)2. For sufficiently small 5, 
we have 

(10) Pi{t\'H\)^\\t\H\)5, i = l,2, and Pl;^{t\'Ht) - X^'^it\'Ht)5^ 

for all t gT- Again following Ventura, Cai and Kass (2005), we may smooth 
the observed- frequency estimates of P^'i^tlTlt) to produce an estimate of 
X^''^{t\T-Lt), and smooth the observed-frequency estimates of Pl{t\7il) to pro- 
duce estimates of X^{t\7il), z = 1, 2. Letting (5 — > in (8), we obtain 

(11) X''\t\'Ht) = m>^Ht\'Hl)X\t\n^) Vte [0,T). 

Consequently, for sufficiently small 6, a conditional test of HQ:^{t) = 1 for 
ah t becomes a test of the null hypothesis Hq : X^'^{t\nt) = X^{t\nl)X^{t\nt) 
for all t or, equivalently, in this conditional case we have the same null 
hypothetical statement as (7). 

In attempting to make equation (9) rigorous, a difficulty arises: for a regu- 
lar marked point process, the function ^ need not be independent of the spik- 
ing history. This would create a fundamental mismatch between the discrete 
data- analytical method and its continuous-time limit. The key to avoiding 
this problem is to enforce the sparsity condition (10). Specifically, the prob- 
abilities PKtlWf) are of order 0{5), while the probabilities Pi'iltlUt) are 
of order 0((5^). This also allows independence models within the marked 
point process framework. Section 6 proposes a class of marked point process 
models indexed by 6 and provides results that validate the heuristics above. 

2.3. Loglinear models. We now reformulate in terms of loglinear models 
the procedures sketched in Sections 2.1 and 2.2 for v = 2 neurons, and then 
indicate the way generalizations proceed when u >3. 

In the marginal case of Section 2.1, it is convenient to define 

pl:^{t)=pht), 
p^:?{t)=pht), 
Pi;iit) = Pi;iit) vter. 

Equation (3) implies that 

(12) log[P'J{t)] = alog[Pl{t)] + 61og[Pf (i)] + ablog[at)] 

for all a, 6 € {0, 1} and t G T and in the continuous-time limit, using (6), we 
write 

(13) log Ai'2 (t) = log A^ (t) + log A^ (t) + log[e(t)] 
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for t S [0,r). The null hypothesis may then be written as 

(14) Fo:log[e(t)]=0 Vte[0,T). 

In the conditional case of Section 2.2, we similarly define 

plfmt)=pimt), 

Pi;i{t\nt) = Pi;i{t\nt) yteT, 

we may rewrite (8) as the loglinear model 

(15) log[P'Jmt)]=alog[PHt\nl)] + blog[P?mt)]+ablog[C{t)] 

for all a, 6 S {0, 1} and t ^T and in the continuous-time limit we rewrite 
(11) in the form 

(16) \og\^^\t\'H,) = log XHm]) + log X^ATii) + log[e(t)] 

for t € [0,T). The null hypothesis may again be written as in (14). 

Rewriting the model in loglinear forms (12), (13), (15) and (16) allows us 
to generalize to z^ > 3 neurons. For example, with the obvious extensions of 
the previous definitions, for i/ = 3 neurons the two-way interaction model in 
the continuous-time marginal case becomes 

^og[Plf!{t)] = log[Pht)] + log[Pf (t)] + log[P!{t)] 
(17) 

-fafelog[C{i,2}(i)] +«clog[C{i,3}(t)] + 6clog[C{2,3}(i)] 

for all a,b,c£ {0, 1} and t ^T, and 

log[Ai'2'3(t)] = log[X\t)] + log[X\t)] + log[A3(t)] 

+ log[e{l,2}(t)] +log[C{l,3}(t)] +log[e{2,3}(i)] 

for all t € (0,T]. The general form of (17) is given by equation (28) in Sec- 
tion 5. In the conditional case, the two-way interaction model becomes 

^og[Pa'bcmt)] = aiog[Pimt)] + biog[pfmt)] + ciogipfmt)] 

(18) 

+ ablog[C{i,2}{t)] + aclog[C{i,3}(t)] + ?)clog[C{2,3}(i)] 

for all a,b,c(z {0,1} and t^T and in continuous time, 

\og[X^^^^\t\Ht)] = \og[X\t\nt)] + \og[X\t\Ht)] + \og[X\t\Ht)] 
+ log[e{i,2}(i)] +log[e{i,3}(i)] +log[^{2,3}(t)] 
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for all t E (0,r]. In either the marginal or conditional case, the null hypoth- 
esis of independence may be written as 

(19) i/o:logK{„}(t)]=0 VtG(0,r],l<i<j<3. 

On the other hand, we could include the additional term afeclog[^M 2,3}(^)] 
and use the null hypothesis of no three-way interaction 

(20) i:^o:log[e{i,2,3}(i)]=0 VtG(0,r]. 

These loglinear models offer a simple and powerful way to study depen- 
dence among neurons when spiking history is taken into account. They have 
an important dimensionality reduction property in that the higher-order 
terms are asymptotically independent of history. In practice, this provides 
a huge advantage: the synchronous spikes are relatively rare; in assessing 
excess synchronous spiking with this model, the data may be pooled over 
different histories, leading to a much larger effective sample size. The gen- 
eral conditional model in equation (34) retains this structure. An additional 
feature of these loglinear models is that time-varying covariates may be 
included without introducing new complications. In Section 3 we use a co- 
variate to characterize the network up states, which are visible in part C of 
Figure 1 , simply by including it in calculating each of the individual-neuron 
conditional intensities X^{t\'Hl) and X'^{t\'H'f) in (16). 

Sometimes, as in the data we analyze here, the synchronous events are 
too sparse to allow estimation of time- varying excess synchrony and we must 
assume it to be constant, C{t) = ( for all t. Thus, for u = 2, the models of (12) 
or (15) take simplified forms in which (^(t) is replaced by the constant (^ and 
we would use different test statistics to test the null hypothesis Hq : ^ = 1 . 
To distinguish the marginal and conditional cases, we replace C(t) by (^h 
in (15) and then also write Hq-.Qh = 1- Moving to continuous time, which 
is simpler computationally, we write ^(t) = ^, estimate ^ and ^h, and test 
Hq:E, = 1 and Hq:S,h = 1- Specifically, we apply the loglinear models (12), 
(13), (15) and (16) in two steps. First, we smooth the respective PSTHs to 
produce smoothed curves A'(t), as in parts A and B of Figure 1. Second, 
ignoring estimation uncertainty and taking A*(i) = A*(i), we estimate the 
constant (. Using the point process representation of joint spiking (justified 
by the results in Sections 5 and 6), we may then write 



log L(0 = - y A(t) dt + J^ log A(ii 



where the sum is over the joint spike times tj and A(t) is replaced by the 
right-hand side of (13), in the marginal case, or (16), in the conditional case. 
It is easy to maximize the likelihood L{^) analytically: setting the left-hand 
side to i{C), in the marginal case we have 



/N 
X\t)X\t)dt + j, 
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where A^ is the number of joint (synchronous) spikes (the number of terms 
in the sum), while in the conditional case we have the analogous formula 

/AT 
x\t\n])\^{t\'Hi)dt + — 

and setting to and solving gives 

(21) i ^ 



and 

N 
(22) in 



which, in both cases, is the ratio of the number of observed joint spikes to 
the number expected under independence. 

We apply (21) and (22) in Section 3. To test Hq-.C = ^ and Hq-.^h = I, 
we use a bootstrap procedure in which we generate spike trains under the 
relevant null- hypothetical model. This is carried out in discrete time, and re- 
quires all 4 cell probabilities P^'^ (t) or P^'^ {t\T^t) at every time t ^T. These 
are easily obtained by subtraction using Pi{t) = X^{t)6, -Pf(t) = X'^{t)S, 
and C = I or, in the conditional case, Plit\nl) = X^itlTiDS, Pl{t\rQ) = 
}?{t\'H'^)5, and Ch = ^H- As we said above, A*(t) = A*(t) is obtained from 
the preliminary step of smoothing the PSTH. Similarly, the conditional in- 
tensities X'^{t\7il) = X^{t\'Hl) are obtained from smooth history-dependent in- 
tensity models such as those discussed in Kass, Ventura and Brown (2005). 
In the analyses reported here we have used fixed-knot splines to describe 
variation across time t. 

In the case of 3 or more neurons the analogous estimates and cell proba- 
bilities must, in general, be obtained by a version of iterative proportional 
fitting. For i^ = 3, to test the null hypothesis (20), we follow the steps leading 
to (21) and (22). Under the assumption of constant Ci23) we obtain 

N 
(23) 623 



/Ai(t)A2(i)A3(i)62(t)ei3(i)63(i)di 
and 

(24) ^123,/f=pi^^|^i^^2(i|^2)^3(t|^3)^^2^^(i)^^3^^(i)^23^^(i)^i- 

In Section 3 we fit (17) and report a bootstrap test of the hypothesis (20) 
using the test statistic ^123 in (23). 



LOGLINEAR MODELS OF SYNCHRONY 15 

3. Data analysis. We applied the methods of Section 2.3 to a subset of 
the data described in Section 1.1 and present the results here. We plan to 
report a more comprehensive analysis elsewhere. 

We took (5 = 5 milliseconds (ms) , which is a commonly-used window width 
in studies of synchronous spiking. Raster plots of spike trains across repeated 
trials from a pair of neurons are shown in Parts A and B of Figure 2, with 
synchronous events indicated by circles. Below each raster plot is a smoothed 
PSTH, that is, the two plots show smoothed estimates A^(t) and A^(t) of 
X^{t) and A^(t) in (6), the units being spikes per second. Smoothing was 
performed by fitting a generalized additive model using cubic splines with 
knots spaced 100 ms apart. Specifically, we applied Poisson regression to the 
count data resulting from pooling the binary spike indicators across trials: 
for each time bin the count was the number of trials on which a spike oc- 
curred. To test Hq under the model in (13), we applied (21) to find log^. 
We then computed a parametric bootstrap standard error of log^ by gen- 
erating pseudo-data from model (13) assuming Hq : log^ = 0. We generated 
1000 such trials, giving 1000 pseudo-data values of log^, and computed the 
standard deviation of those values as a standard error, to obtain an observed 
z-ratio test statistic of 3.03 {p = 0.0012 according to asymptotic normality). 

The highly significant z ratio shows that there is excess sychronous spiking 
beyond what would be expected from the varying firing rates of the two 
neurons under independence. However, it does not address the source of the 
excess synchronous spiking. The excess synchronous spiking could depend on 
the stimulus or, alternatively, it might be due to the slow waves of population 
activity evident in part (C) of Figure 1, the time of which vary from trial 
to trial and therefore do not depend on the stimulus. To examine the latter 
possibility, we applied a within-trial loglinear model as in (16) except that we 
incorporated into the history effect not only the history of each neuron but 
also a covariate representing the population effect. Specifically, for neuron i 
{i = 1,2) we used the same generalized additive model as before, but with 
two additional variables. The first was a variable that, for each time bin, was 
equal to the number of neuron i spikes that had occurred in the previous 
100 ms. The second was a variable that, for each time bin, was equal to 
the number of spikes that occurred in the previous 100 ms across the whole 
population of neurons, other than neurons 1 and 2. We thereby obtained 
fitted estimates X^itlV.}) and \^{t\'Hl) of \^{t\'H\) and \^{t\H'l). Note that 
the fits for the independence model, defined by applying (7) to (11), now 
vary from trial to trial due to the history effects. Applying (22), we found 
log^H-, and then again computed a bootstrap standard error of log.^^ by 
creating 1000 trials of pseudo-data, giving log^ji^ = 0.06 ± 0.15, for a z-ratio 
of 0.39, which is clearly not significant. 

Raster plots for a different pair of neurons are shown in parts (E) and (F) 
of Figure 2. The same procedures were applied to this pair. Here, the z-ratio 
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for testing Hq under the marginal model was 3.77 {p < 0.0001), while that for 
testing Hq under the conditional model remained highly significant at 3.57 
{p = 0.0002) with log^H = 0.82 it 0.23. In other words, using the loglinear 
model methodology, we have discovered two pairs of VI neurons with quite 
different behavior. For the first pair, synchrony can be explained entirely 
by network effects, while for the second pair it can not; this suggests that, 
instead, for the second pair, some of the excess synchrony may be stimulus- 
related. 

We also compared the marginal and conditional models (13) and (16) 
using ROC curves. Specifically, for the binary joint spiking data we used 
each model to predict a spike whenever the intensity was larger than a given 
constant: for the marginal case whenever logA^'^(t) > Cmarginah and for the 
conditional case whenever logX^'^{t\T-Lt) > Cconditionai ■ The choice of constants 
Cmarginai and Cconditionai reflect trade-offs between false positive and true posi- 
tive rates (analogous to type I error and power) and as we vary the constants, 
the plot of true vs. false positive rates forms the ROC curve. To determine 
the true and false positive rates, we performed ten-fold cross-validation, re- 
peatedly fitting from 90% of the trials and predicting from the remaining 
10% of the trials. The two resulting ROC curves are shown in part D of Fig- 
ure 2, labeled as "no history" and "history," respectively. To be clear, in the 
two cases we included the terms corresponding, respectively, to ^ and S,h, 
and in the history case we included both the auto-history and the network 
history variables specified above. The ROC curve for the conditional model 
strongly dominates that for the marginal model, indicating far better pre- 
dictive power. In part C of Figure 2 we display the true positive joint spike 
predictions when the false-positive rate was held at 10%. These correctly- 
predicted joint spikes may be compared to the complete set displayed in 
parts A and B of the figure. The top display in part C, labeled "no history," 
shows that only a few joint spikes were correctly predicted by the marginal 
model, while the large majority were correctly predicted by the conditional 
model. Furthermore, the correctly predicted joint spikes are spread fairly 
evenly across time. In contrast, the ROC curves for the second pair of neu- 
rons, shown in part (G) of Figure 2, are close to each other: inclusion of 
the history effects (which were statistically significant) did not greatly im- 
prove predictive power. In (G), the correctly predicted synchronous spikes 
are clustered in time, with the main cluster occurring near a peak in the 
individual-neuron firing-rate functions shown in the two smoothed PSTHs 
in parts (E) and (F). 

Taking all of the results together, our analysis suggests that the first pair 
of neurons produced excess synchronous spikes solely in conjunction with 
network effects, which are unrelated to the stimulus, while for the second 
pair of neurons some of the excess synchronous spikes occurred separately 
from the network activity and were, instead, stimulus-related. 
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Fig. 3. Plots of synchronous spiking events for 3 neurons. Each of the three plots displays 
all joint spikes (as circles) for a particular pair of neurons. The dark circles in each plot 
indicate the 3-way joint spikes. 



We also tried to assess whether 2-way interactions were sufficient to ex- 
plain observed 3-way events by fitting the no-3-way interaction model given 
by (17), and then testing the null hypothesis in (20). We did this for a par- 
ticular set of 3 neurons, whose joint spikes are displayed in Figure 3. The 
method is analogous to that carried out above for pairs of neurons, in the 
sense that the test statistic was ^123 given by (23) and a parametric bootstrap 
procedure, based on the fit of (17), was used to compute an approximate 
p-value. Fitting of (17) required an iterative proportional fitting procedure, 
which we will describe in detail elsewhere. We obtained p = 0.16, indicating 
no significant 3-way interaction. In other words, for these three neurons, 2- 
way excess joint spiking appears able to explain the occurrence of the 3-way 
joint spikes. However, as may be seen in Figure 3, there are very few 3-way 
spikes in the data. We mention this issue again in our discussion. 

4. A marked point process framework. In this section a class of marked 
point processes for modeling neural spike trains is briefly surveyed. These 
models take into account the possibility of two or more neurons firing in syn- 
chrony (i.e., at the same time). Consider an ensemble of i' neurons labeled 1 
to u. For T > 0, let Nt denote the total number of spikes produced by this 
ensemble on the time interval [0, T) and let < si < • • • < sn^ < T denote 
the specific spike times. For each j = 1, . . . , Nt, we write (sj, {ii, . . . ,ik)) to 
denote the event that a spike was fired (synchronously) at time Sj by (and 
only by) the ii, . . ., i^ neurons. We observe that 



(25) 



'Ht = {(si, Ki), . . . , {snt,i^Nt) : Kj G /C, j = 1, . . . , Nt} 



18 



R. E. KASS, R. C. KELLY AND W.-L. LOH 



forms a marked point process on the interval [0, T) with /C as the mark space 
satisfying 

/CC {{ii,...,ik):l<ii <■■■ <ik<i^,k = l,...,i^}. 

We follow Daley and Vere-Jones (2002), page 249, and define the conditional 
intensity function of T-Lt as 

\{t,K\'Ht) 

= A(t, k|{(si, Ki), . . . , {sNt,KNt)}) 

(hl{t)fl{K\t) VO<t<Sl, 

hi{t\{si,Kl), ..., (Si_l, Ki_l))/j(K|t; (Sl, Kl), . . . , (Sj_l, Kj_i)) 

Vsi_i <t<Si,i = 2,..., Nt, 

hNr+liAi^l^l'^l)^ • • • 1 {sNt,I^NT))fNT+l{f^\'l^'i isi,Ki),..., (sATy, KaTj,)) 

Vsatj, <t <T, 

where /ii(-) is the hazard function for the location of the first spike si, 
h2{-\isi,Ki)) the hazard function for the location of the second spike S2 
conditioned by (si, ki), and so on, while fi{-\t) is the conditional probability 
mass function of ki given si = t, and so on. It is also convenient to write 
A(t, k|0) = X{t,K.\7it) for all t < si. The following proposition and its proof 
can be found in Daley and Vere-Jones (2002), page 251. 



by 



Proposition 1. Let Ht be as in (25). Then the density ofTir is given 



PxCHt) =Px{{{si,ki),...,{snt,i^Nt)}) 



Nt 



JjA(Sj,Kj|?^sJ 



.i=l 



exp 



V / x{t,K\nt)dt 



5. Theoretical results: Marginal methods. In this section we (i) provide 
a justification of the limiting statements in (6) and (ii) generalize to higher- 
order models. We also note that lagged dependence can be accommodated 
within our framework, treating the case v = 2. 

5.1. Regular marked point process and loglinear modeling. In this sub- 
section we prove that the heuristic arguments of Section 2 for marginal 
methods hold under mild conditions. Consider y>l neurons labeled 1 to z^. 
For T > 0, let Nt denote the total number of spikes produced by these v 
neurons on the time interval [0, T) and let < si < • • • < snt < T denote the 
specific spike times. For each j = 1, . . . , Nt, we write {sj, (ij)) to represent 
the event that a spike was fired at time Sj by neuron ij where ij € {1, . . . , z^}. 
We observe from Section 4 that 

T-l-T = {(si, {ii)), ■■■, {snt, iiNr))} 
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forms a marked point process on the interval [0, T) with mark space /C = 
{(1), . . . , {i')}. Following the notation of Section 4, let X{t, {i)\'Ht) denote the 
conditional intensity function of the point process T-Lt- We assume that the 
following two conditions hold: 

Condition (I). There exists a strictly positive refractory period for each 
neuron in that there exists a constant ^ > such that A(t, (i) \T-Lt) = if there 
exists some (s^(i)) G Tit such that t — s<9,i&{l,... ,i^}. 

Condition (II). For each k £ {0,. . . ,2\T/9~\ —1} and i,ii, . . . ,ifc € {1, . . . , z^}, 
the conditional intensity function A(t, {i)\{{si, (ii)), . . . , (s^, (ik))}) is a con- 
tinuously differentiable function in (si, . . . , Sk, t) over the simplex < si < 
■■■<Sk<t<T. 

li 6 <9, then condition (I) implies that there is at most 1 spike from each 
neuron in a bin of width 6. Conditions (I) and (II) also imply that the marked 
point process is regular in that (exactly) synchronous spikes occur only with 
probability 0. Theorem 1 below gives the limiting relationship between the 
bin probabilities of the induced discrete-time process and the conditional 
intensities of the underlying continuous-time marked point process. 

Theorem 1. Suppose that conditions (I) and (II) hold, 1 <ii < ■ ■ ■ < 
ik^f and \ <k <v. Then 

\\u.5-^Pi'-t{t^) 

k 

^ E i?n[A(t,(i,)i{(t,(ii))} 



k\ 

iivjfc:{iiv,ife}={n,.--,«fc} '=1 



U---U{(t,(jfe_i))}U7^i)], 



where tm = rnS — )■ t as 5 — >■ and \{t, {i2)\{(t, (ii))} U Ht) = lim.t*^t- A(i, 
(^2)|{(i*5 (^i))} U'Ht), etc. Here the expectation is taken with respect to T-Lf 

Theorem 1 validates the heuristics stated in (6) where i/ = 2, 

\\t) = E[\{t,{i)\nt)i 
X''\t) = 1 Y. ^[^(*' (^2)l{(t, (n))} u nt)x{t, {^^)\nt)]. 

Next we construct the discrete-time loglinear model induced by the above 
marked point process. First define recursively for t^ = m5, 

C{h}itm) = 6-^Pt (tm) Vii = 1, . . . , zy, 

6-^Pl'f^t^) 
C{n,i2}(im) = 7 7-^7 7-^ Vl<^l<^2<^/, 

^{ii}{tm)(,{i2}{tm) 
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(26) 



C{ii,...,«s,}(*m) 



nE:c{Ji,...,ifc}CH(im) 

VI < ii < • • • < ifc < !/, 2 < A; < I/. 
We further define 

(27) e{i,,„„,,}(t) = limC{i,....,^,}(tm) Vl<ii<---<ifc<i/,l</c<zy, 

where hm5_^o tm — ^ t, whenever the expression on the right-hand side of (27) 
is weh defined. The following is an immediate corollary of Theorem 1. 

Corollary 1. Let S,{i^}{t) and (,^i-^i^y{t) be as in (27). Then with the 
notation and assumptions of Theorem 1, we have 

ei,,l(t)=i?[A(t,(n)|?^t)], 

SC{ii,...,ife} 

k 

X E E\{[\{t,{hmt,{h))} 

{3l,-,3k}={n,...,ik} ^=1 

U---U{(t,(jfc_i))}U?^i)], 
whenever the right-hand sides are well defined. 

It is convenient to define PQ'"''Q{tm) = 1- For ai, . . . ,0,^ G {0, 1} and not 

all 0, define Pal','.'.'!^a^{tm) = Pi^'"'i'' (tm) where z G {zi, . . . ,ik} if and only if 
Oi = 1. Then using the notation of (26), the corresponding loglinear model 
induced by the above marked point process is 

loglP^-'^aAt^)] 

(28) =log[P;Xi''(*-)] 

= Y,(^^^Og[Pl{tm)] + Y. (n«j) l°g[CH(im)] 

i=l HC{1,...,!/}:|H|>2 S'eS ^ 

for all ai, . . . ,au € {0, 1}. Under conditions (I) and (II), Corollary 1 shows 
that ^sit) = lim^-s-o CB.{tm) is continuously differentiable. This gives an asymp- 
totic justification for smoothing the estimates of Cs, ^ ^ {1, • • • , i^}- 
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5.2. Case of v = 2 neurons with lag h. This subsection considers the 
lag h case for two neurons labeled 1 and 2. Let h,m be integers such that 
Q <m<m + h< T5~^ — 1. As in (1), we write 

Pa:bitm,tm+h) = P[X\t^^h) = b, X\tm) = a] Va, 6 € {0, 1}, 

where tm = m5 and tm+h = {"m- + h)6. Analogous to Theorem 1, we have the 
following results for the lag case. 

Theorem 2. Suppose conditions (I) and (II) hold. Then 

P,:,{tm,tm+H) ^ ^^^^^ ^ ^ ^2)|{(t, (1))} U nt+r)Kt, (1)1^*)], 
5->-0 0^ 

where t^+h -^ t + r and t^ —^t as 6 —^0 for some constant r > 0. Here the 
expectation is taken with respect to Tit+T (and hence also lit)- 

Corollary 2. Let C,{tm,tm+h) be defined as in (5). Then with the no- 
tation and assumptions of Theorem 2, we have 

\\mC,{tm,tm+h) 
(5— !>0 

^^^^ E[x{t + T,{2)\{{t,{i))}unt+r)x{t,{i)\nt)] vn<f^T- 

E[\{t + T,{2)\'Ht+r)]E[\{t,{l)\nt)] 
whenever the right-hand side is well defined. 

We observe from conditions (I) and (II) that the right-hand side of (29) 
is continuously differentiable in t. Again this provides an asymptotic justifi- 
cation for smoothing the estimate of C,{t, t -|- r), with respect to i, when 5 is 
small. 

6. Theoretical results: Conditional methods. This section is analogous 
to Section 5, but treats the conditional case. We (i) provide a justification 
of the limiting statements in (9) and (ii) generalize to higher-order models. 
We again also note that lagged dependence can be accommodated within 
our framework, treating the case v = 2. 

6.1. Synchrony and loglinear modeling. This subsection considers z^ > 1 
neurons labeled 1 to z^. We model the spike trains generated by these neurons 
on [0, r) by a marked point process "Ht with mark space 

/C = {(zi,...,Zfc):l <ii<--- <ik<u,k = l,...,i^}. 

Here, for example, the mark (ii) denotes the event that neuron ii (and only 
this neuron) spikes, (^1,^2) denotes the event that neuron ii and neuron 12 
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(and only these two neurons) spike in synchrony (i.e., at the same time), and 
the mark (1, . . . , z^) denotes the event that all v neurons spike in synchrony. 
Let iVj denote the total number of spikes produced by these neurons on 
[0,i), < t < T, and let < si < • • • < satj, < T denote the specific spike 
times. For each j = 1, . . . , A't^, let kj G /C be the mark associated with Sy 
Then Ht can be expressed as 

(30) •HT = {(si,Ki),...,(s7Vy,KAry)}. 

Given "Ht, we write 
"^t = {•^ • (S)^) ^ ^t fo^ some K = (ii, . . . ,ik) such that i G {ii, . . . ,ifc}} 

yi = l, . . . ,1/. 

Til denotes the spiking history of neuron i on [0, t). The conditional intensity 
function A(t, K.\7it), t G [0,T) and k G /C, of the marked point process "Ht is 
defined to be 

xit,{i)\nt) = x'it\7{i) vtG[o,r), 

(31) 

A: 

A(t, (zi, . . .,ik)\nt) = <5'^-S{n,...,*.}(i) n ^''(t\^') V* e [o,r), 

where 5 > is a constant, 7{ji,...,ij.}(i)'s are functions depending only on t 
and the A*(i|?^^)'s are conditional intensity functions depending only on the 
spiking history of neuron i. We take Ju\{t) to be identically equal to 1. 

From (31), we note that the above marked point process model is not 
a single marked point process but rather a family of marked point processes 
indexed by 5. In the sequel, we let 5 — t- 0. We further assume that the fol- 
lowing two conditions hold: 

Condition (III) . There exists a strictly positive refractory period for each 
neuron in that there is a constant ^ > such that, for i = 1, . . . , z/ and t G 
[0,T), 

A*(i|'Hj) = 0, if there exists some s eV-l such that t — s <9. 

Condition (IV). For each /c G {0, . . . , \T/e] - 1} and i G {1, ... , u}, X'{t\{si, 
...,Sfc}) is a continuously differentiable function in (si, . . . ,Sfc,t) over the 
simplex < si < • ■ ■ < s^ < i < T. 

Following Section 2.2, we divide the time interval [0,T) into bins of 
width 6. For simplicity, we assume that T is a multiple of 6. Let t^ = md 
and X'{tm), m = 0, . . . ,T6-^ - 1, be as in Section 1. If X^{ti) = 1 for ah 
I € {li, . . . , Ik}, and X'^iti) = otherwise, for some subset <li < ■ ■ ■ <lk < 
m — 1, we write 

(32) 7il^ = {tij^,...,tii^}, 7it^ = {'Ht^,...,7it^). 
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It should be observed that although the above definitions of T-Li and T-Li 
differ from those given in Section 1.2, they are equivalent. We note that 
the conditional intensity functions A(t, (ii, . . . ,«fc)|^t) in (31) depend on the 
bin width 5. This is necessary in order to preserve the natural hierarchical 
sparsity conditions given by 

sup P[X^l(im) = l,...,^*H*m) = l]=0(5'=), 

m=0,...,T<5-i-l 

as (5 —7- for all 1 < ii < • • • < ^fc ^ ^) 1 < A; < z/. 

Theorem 3. Consider the marked point process Tix as in (30) with 
conditional intensity function satisfying (31). Then under conditions (III) 
and (IV), we have 

PKtm\ntj = 6x\tm\niJ + o{s''), 

2 
Plt(*rn]ntJ = <52[1 +7{n..}(tn^)] ll^>^'^ {tm\V^J] + 0{6'), 

and in general, 

k J=l Hi,...,Hj:all disjoint and noncmpty,UHj={ii,...,jfe} i=l J 



k 
X 



ll[x''{tm\ni^J] + o{6'+'] 



for sufficiently small 6 where T-Ltm o.nd Til^ are defined by (32). 

The following is an immediate corollary of Theorem 3. It gives an asymp- 
totic justification for equation (8) in Section 2.2. 

Corollary 3. With the notation and assumptions of Theorem 3, we 
have for v = 2, 

Pli{tm\ntJ = C{tm)Pl{tm\nl)P?{t„,\nl) + 0(6') 

for sufficiently small 6 > uniformly over Hl^, m = 0, . . . ,T6~^ — 1 where 

C(tm) = l +7(1,2} (im). 

We now use Theorem 3 to construct a loglinear model (for the above spike 
train data) whose higher-order coefficients are asymptotically independent 
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of past spiking history. First define recursively 

C{h,i2]{tm\HtJ = ^ ^ 'il<li<l2<U, 

^{h}\tm\rtt^}(^{i^{tm\n.t,^) 

(33) 

C{^,,...,i^}{tm\ntJ = ^ . ., lay N Vf < ii < . . . < ifc < I/. 

It follows from Theorem 3 and (33) that for sufficiently small 5, 
C{n,i2}(*-I^t™) = 1 + l{h,i2}itm) + 0{5) VI < ii < i2 < y, 



C{Ji,...,ifc}(im|^tm) 

2-^j = l Z_/Hi,...,Hj:all disjoint and nonempty,USj ={ii ,. . .,ifc} 1 li = l '"i v^ 

nE:c{ii,...,ifc}:| = |>2CH(im|^t™^) 



+ 0(5) 



whenever 1 < ii < • • • < ifc ^ i^ and k>2, assuming that terms on the right- 
hand side are well defined. The practical importance of these results lies in 
the fact that the coefficients C{ji,...,ife}(*m|^tm) with k>2 are asymptotically 
(as (5 — > 0) independent of T^t^ , the spiking history of the neurons. It is 
convenient to define -Pq' 'o'(*m|'^t„) = 1. For ai, . . . , a,y G {0, 1} and not all 0, 

define -Par;.'.'.ra.(im|^t™) = -Pi'.',',',i**(^m|'^t™) where i € {ii, . . . ,ifc} if and only 
if Oj = 1. Then the induced loglinear model is 

^og[p^l■■■::a^t„^\'HtJ] 

(34) 

= ^aiiog[Piitm\ntj] + E (n °^) ^°smtm\ntj] 

i=l =C{1,...,i/}:|H|>2S'GH ^ 

for all ai, . . . , Oj/ G {0, 1} where the second term on the right-hand side of (34) 
is asymptotically (as (^ — )• 0) independent of the spiking history Tltm ■ 

6.2. u = 2 neurons with time- delayed synchrony. This subsection consid- 
ers z^ = 2 neurons labeled 1, 2 and let r > be a constant denoting the spike 
lag. We model the spike train generated by the two neurons on [0,T) by 
a marked point process "Hj- as in (25) with mark space /C = {(1), (2), (1, 2)}. 
The marks (1), (2) are interpreted as before as isolated (i.e., nonsynchronous) 



LOGLINEAR MODELS OF SYNCHRONY 25 

spikes. However, now (1,2) is interpreted to be neuron 1 spiking first and 
then neuron 2 spiking second after a delay of r time units. The mark (1, 2) is 
used to model a precise time-delayed synchronous spiking of lag r between 
the 2 neurons. 

Let Nt denote the number of times the three marks occur on [0,T) and 
si < • • • < snt be the specific spike times. For each j = 1, . . . , Nt, let Kj G /C 
be the mark associated with Sj. Then T-Lt can be decomposed into {Ti.} ^T-LIj^^-) , 
where 

n\ = {s : (s, k) G Tit where k G {(1), (1, 2)}}, 
(35) 

n^t+^ = {s:{s,K)ent+rwheieKe{{2),{l,2)}}. 

To be definite, {s,k) = {s, (1,2)) means neuron 1 spikes at time s and neu- 
ron 2 spikes at time s + t. The conditional intensity function A(t,K|?^4), 
t G [0,T) and k G /C, of the marked point process Ht is defined to be 

x{t,{i)\nt) = x\t\nD vi = i,2, tG[o,T), 

(36) 

x{t,{i,2)\nt) = 6jit,t + T)\\t\'Hl)x\t + T\n'^^^) vtG[o,r), 

where 5 > is a constant, 7(t, t + r) is a continuously differentiable function 
in t on the interval < t < T — r, and X^{t\'Hl), X^ {t + T\T-lf_^^) are conditional 
intensity functions depending only on the spiking history of neuron 1 up to 
time t and on the spiking history of neuron 2 up to time t + T, respectively. 
As in Section 1, we divide the time interval [0, T) into bins of width 5 > 0. 
Let Htm — (^tm'^tm+h) ^^ ^^ ™ ("^^^ where, for simplicity, we assume that 5 
is chosen such that m, h are integers satisfying Q<m<m + h< T5~^ — 1 
and tm+h = tm + T. Recall that, by definition, 

Pl{tm\niJ = P{x\tm) = i\nij, 

Pl^^{tra,tm+h\'HtJ = P{X\tm) = l,X\tra+h) = l|^t„). 

Theorem 4. Consider the marked point process T-Lt ols in (35) with 
conditional intensity function satisfying (36). Let m, h be integers satisfying 
<m <m-\- h < T6~^ — 1 and tm+h =tm + T. Then under conditions (III) 
and (IV), we have 

Pl'^{tm, tm+h\nt^ ) = S^ Mtm., tm+h) + 1] A' (^m |^L ) A' {tm+h\n}^^, ) + 0{5^) 

for sufficiently small 5. 

The practical significance of Theorem 4 is that 'y{tm,tm+h) does not de- 
pend on the spiking history of the 2 neurons. This implies that a statistic 
based on 7 can be constructed to test the null hypothesis Hq that there is 
no time-delayed spiking synchrony at lag r between the 2 neurons. 
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7. Discussion. We have described an approach to assessing spike train 
synchrony using loglinear models for multiple binary time series. We tried 
to motivate the application of loglinear modeling technology in Section 1, 
emphasizing two features of individual neural response: stimulus-induced 
nonstationarity that remains time-locked across trials, and within-trial ef- 
fects that are history-dependent, with timing that varies across trials. These 
were incorporated into the models by including for individual neurons both 
time- varying marginal effects, which stay the same across trials, and history- 
dependent terms; interaction terms were treated separately. In Section 3 we 
presented results for two pairs of neurons. For both pairs there was evidence 
of excess synchronous spiking beyond that explained by stimulus-induced 
changes in individual-neuron firing rates. In one pair, network activity, rep- 
resented as history dependence, was sufficient to account for excess syn- 
chronous spiking, but the other pair displayed excess synchronous spiking 
that remained highly statistically significant even after network effects were 
incorporated, indicating stimulus-related synchrony. Our theoretical results 
provided a continuous-time point process foundation for the methods, justi- 
fying both our use of smoothing and our derivation of the excess-synchrony 
estimators ( and Ch- 

Assessment of synchrony via continuous-time loglinear models is closely 
related to the unitary-event analysis of Griin, Diesmann and Aertsen (2002a, 
2002b). Unitary event analysis assumes each neuron follows a locally-statio- 
nary Poisson process, which has been shown to be somewhat conservative 
in the sense of providing inflated p-values in the presence of non-Poisson 
history dependence. Its main purpose is to identify stimulus-locked excess 
synchrony. Because the loglinear models could be viewed as generalizations 
of locally-stationary Poisson models, they could extend unitary-event analy- 
sis to cases in which it seems desirable to account more explicitly for stimulus 
and history effects. This is a topic for future research. 

We also provided an example of testing for 3-way interaction. The results 
we gave in Section 3 for a particular triple of neurons indicated no evidence 
of excess 3-way joint spiking above that explained by 2- way joint spiking. 
A systematic finding along these lines, examining large numbers of neurons, 
would be consistent with findings of Schneidman et al. (2006). However, as 
may be seen from Figure 3, 3-way joint spikes are very sparse. A careful 
study of the power to detect 3-way joint spiking in contexts like the one 
considered here could be quite helpful. We plan to carry out such a study 
and report it elsewhere. 

We have restricted history effects to individual neurons by assuming, first, 
that each neuron's history excludes past spiking of the other neurons un- 
der consideration and, second, that the interaction effects are independent 
of history. This greatly simplifies the modeling and avoids confounding the 
interaction effects with cross-neuron effects. While it would be possible, in 
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principle (by modifying the hierarchical sparsity condition) , to allow history- 
dependence within interaction effects, we see no practical benefit of doing 
so. With the two additional, highly plausible assumptions used here, we get 
both tractable discrete-time methods and a sense in which the methods may 
be understood in continuous time. A key element of our formalism is the re- 
quirement of hierarchical sparsity, as in the special case of equation (10) and 
more generally in Section 6.1 (preceding Theorem 3). This corresponds to the 
practical reality that two-way synchronous spikes are rare, as in Figure 2, 
and three-way spikes are even rarer, as in Figure 3. Some form of spar- 
sity seems to us essential. [After our article was accepted we became aware 
that Solo (2007) had attempted to develop likelihoods for point processes 
having synchronous events, but because his approach does not incorporate 
sparsity we have been unable to understand how it could be used in the 
kind of applications we have described here.] It is somewhat inelegant to 
have a sequence of marked processes (indexed by 5), but this appears to be 
the best that can be achieved by starting with very natural discrete-time 
loglinear models. An alternative would be to use more standard point pro- 
cess models with short time-scale cross-neuron effects. Presumably, similar 
results could be obtained, but the relationship between these different ap- 
proaches is also a subject for future research. A quite different technology 
involves permutation-style assessment via "dithering" or "jittering" of indi- 
vidual spike times [cf. Geman et al. (2010), Griin (2009)]. Synchrony is one 
of the deep topics in computational neuroscience and its statistical identi- 
fication is subtle for many reasons, including inaccuracies in reconstruction 
of spike timing from the complicated mixture of neural signals picked up 
by the recording electrodes [e.g., Harris et al. (2000), Ventura (2009)]. It 
is likely that multiple approaches will be needed to grapple with varying 
neurophsyiological circumstances . 

APPENDIX 



Proof of Theorem 1. For simplicity, we shall consider only the case 
for 1/ = 2. The proof for other values of i^ is similar though more tedious. We 
observe from Proposition 1 that 

-^1,0 V^m) 



2[T/e] 
6 



1 ' ' . /Tm /Tm. /-tm+l 

E E / ■•■/ / 

fc=o ii,...,ifee{i,2}-^° Jsk-iJt„, 

(37) A(sfc+i,(l)|{(si,(zi)),...,(sfc,(ifc))}) 
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■ k 

llx{si,{iims,,{h)),...,{si.i,{ii.i))}) 



.1=1 



Condition (I) implies that the summations ^^ ^^ ^ g m 2} in (37) contain 
a finite number of summands. Hence, letting (5 — )■ 0, the right-hand side 
of (37) equals 

I ' I 1 ftm ftm ftm + l 

E E 1™,^/ ■■■ / A(s,+i,(i)i{(si,(n)),...,(..,(^.))}) 



(38) X 



nA(s^(^oi{(^i,(n)),...,(5^_i,(^^_i))}) 



i=l 



X e-^Ulo-^' 'i-''i^)\^-)'^-' ds,^,ds,--- ds^. 

Using Condition (II) and the Taylor expansion, we have 

A(sfe+i, (l)|{(si, (ii)), . . . , (sfc, (ife))}) 

= A(t^, (l)|{(si, (ii)), . . . , {sk, {ik))}) + 0{6), 

uniformly over tm < s^+i < tm+i,0 < si < ■ ■ ■ < Sk <tm. Consequently, (38) 
equals 

E E 1™ / •••/ A(i™,(l)|{(si,(n)),...,(Sfc,(Zfc))}) 

A: 



X 



nA(s/,(i/)i{(5i,(ii)),...,(sz_i,(iz_i))}) 

= i?[A(t,(l)|?^i], 
since im+i — ^m = (^ and lim5_j.otm = t- Using a similar argument, we have 



-^1,1 (^mj 

6^~ 



2^ A2 2^ 



§2 Z^ 1^ 

1<^1^£2<2 fc=0 ii,...,ifee{l,2} 



tm rtm rtm + l rtm + 1 

J Si^_^ Jtm -Js^ + i 
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A(sfc+2, (2)|{(si, (ii)), . . . , (sk, (ifc)), (sfe+i, (1))}) 
X A(sfc+i, (l)|{(si, (ii)), . . . , (sfc, (ik))}) 

' k 

X n '^(^'' (^oi{(si> (n)), ■ • ■ , (si-i, («i-i))}) 



.1=1 



^2 r'm+l 



X e-S.^-/o™^^AKO-)|K.)..^,^^^^,^_^^^,^, ... ^^^ 



1 



{i?[A(i,(2)|{(t,(l))}U?^i)A(t,(l)|?^t)] 



+ i?[A(i,(l)|{(t,(2))}U?^i)A(t,(2)|?^0]} 
and <5-ipi(t„) = 6~'[Pl;^{tra) + Pi;'(t™)] ^ E[X{t, {l)\nt] as 5 ^ 0. D 

Proof of Theorem 2. For simplicity, define 
for all a, 6, c, (i, S {0, 1}. We observe from Proposition 1 that 

r-2pl,2;l,2/, , N 

1,0;0;1 V ™' "1+^/ 



^ 2\T/&^2\T/e^ 

j=0 fe=0 ii,...,jj,jj+2,---,ij+fc+i€{l,2} 



E 



C7?i /"tm /"Cm + l r^m + h r^m + h /'^7?i + /i4-l 



_i Jt 



Sj + l ■J Sj^i^ Jtm + h 



Ksj+k+2, (2)|{(S1, (ii)), . . . , (Sj, (ij)), (Sj+1, (1)), (Sj+2, (ii+2)), 

...,(Sj+fc+l,(Zj+fc+l))}) 

"i+fc+1 
X n Ksh{k)\{{si,{h)),...,{sj,{ij)), 

{sj+i, (1)), {sj+2, {ij+2)),- ■ ■ , (s;_i, (i;-i))}) 
X A(sj+i, (l)|{(si, (ii)), . . . , {sj, (ij))}) 

X n '^(^'' (^OI{(si, (n)), ■ • ■ , {si-i, (ii^i))}) 
.1=1 

X e- ELi /o'"+' a(«,,(Q|w.) d«, ^^^.^^^^ ds,+k+i ■ ■ ■ dsj+2 dsj+i ds, ■ 

E[\{t + T, (2)|{(t, (1))} U Ut)\{t, {i)\nt)\ 



■ ■ dsi 
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as 5 -> 0. Theorem 2 follows since Pi'iitm,tm+h) ^ Pi'o'o'i{tm,tm+h) as 
5^ 0. D 

Proof of Theorem 3. For simplicity, we only consider the case z/ = 2. 
Let ni^={ti^-^,...,ti^^J and \Ji=i{h,i, ■ ■ ■ ,k,k,} = {h, ■ ■ ■ Jk}, where 0< 
/i < • • • < 4 < m- 1. We observe from (31) that P/(im|^t„) = ^1^,0 (im,|^t„) + 
0{6^), and 

Pl:oitm\ntjp{ntj i + o{5) f''i+^ /■*'.+! ['-+' 



A^ (sfc+i I {sh.i , • • • , Sii,fe^ }) dsk+i dF{si^ ,---,siJ 
= >^\tm\H]jP{fitJ + 0{6)P{HtJ 

for sufficiently small 5 where F denotes the distribution function of s^^ , . . . , s;^, . 
Here si^ € [ti^,ti^^i), s/^ ^ G [t/^ j,t/^ i+i); etc. This proves the first statement 
of Theorem 3. Next we observe that 

plfitmintjPi-Htj 



52 



^^ '^ J\ Jt. 



h{sk+i)>^^{sk+i\{sh,^ , . . . , SZi,,.^ }) 



X A^(sfc+i |{sz2,i , . . . , szj ^^2 }) c?Sfc+i dF{si^ , • ■ • , s« J 
l + 0(,5) Z-^+i r\+i rt..+i rUr.+i 



<^2 



<!, Jt 



X A^(sfc+i|{si2 1 , . . . , 5(2 j^^ }) (isfc+2 dsk+i dF{si^,. ■■,81,^) 

= [i + 7{trn)]x\trn\nlJx\t^\nl)Pintj + oi6)Pintj 

for sufficiently small 6. This proves Theorem 3. D 

Proof of Theorem 4. We observe that Tit^ -> Ht,,^ is a many-to-one 
mapping and from (36) that 

(39) =6S{t^,tm+h)xHtm\nljx\tm.+h\nl^j 

for sufficiently small 5. We further observe that 

PhUnlj = P{x\tm) = i.x\t^+t,) = o\nlj + 0(6^) 
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= 6X\tm\nl) + 0{6^), 
(40) 

P2(t^+,|7^2^^J =P(Xi(t^) =0,x2(t„,+,) = l\nl^^J + 0i6^) 

Thus, it follows from (39) and (40) that 

P{X\t„^) = l,X\tm+h) = l\'HtJ 

= E[P{X\Un) = l.X\tm+h) = mtj\ntj 

= 5'[7(im, Wh) + l]A'(im|^L)A'(Whl^L+J + 0('5') 

for sufficiently small 5. This proves Theorem 4. D 
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